ON A HIGH-FIDELITY HIERARCHICAL APPROACH 
TO BUCKLING LOAD CALUCATIONS 


JIIAA AAn-4 -4 AAA 

MIMM-iOlU 1- IO^ 


Johann Arbocz* 

Delft University of Technology, The Netherlands 
and 

James H. Starnes** and Michael P. Nemeth*** 
NASA Langley Research Center, Hampton, VA 23681-001 


ABSTRACT 


INTRODUCTION 


As a step towards developing a new design 
philosophy, one that moves away from the 
traditional empirical approach used today in 
design towards a science-based design 
technology approach, a recent test series of 5 
composite shells carried out by Waters [1] at 
NASA Langley Research Center is used. It is 
shown how the hierarchical approach to 
buckling load calculations proposed by Arbocz 
et al [2] can be used to perform an approach 
often called “high fidelity analysis”, where the 
uncertainties involved in a design are 
simulated by refined and accurate numerical 
methods. The Delft Interactive Shell DEsign 
COde (short, DISDECO) is employed for this 
hierarchical analysis to provide an accurate 
prediction of the critical buckling load of the 
given shell structure. This value is used later 
as a reference to establish the accuracy of the 
Level-3 buckling load predictions. As a final 
step in the hierarchical analysis approach, the 
critical buckling load and the estimated 
imperfection sensitivity of the shell are verified 
by conducting an analysis using a sufficiently 
refined finite element model with one of the 
current generation two-dimensional shell 
analysis codes with the advanced capabilities 
needed to represent both geometric and 
material nonlinearities. 
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It is generally agreed that, in order to make 
the development of the Advanced Space 
T ransportation System a success and to 
achieve the very ambitious performance goals 
(like every generation of vehicles lOx safer and 
lOx cheaper than the previous one), one must 
make full and efficient use of the technical 
expertise accumulated in the past 50 years or 
so, and combine it with the tremendous 
computational power now available. It is 
obvious that with the strict weight constraints 
used in space applications these performance 
goals can only be achieved with an approach 
often called “high fidelity analysis", where the 
uncertainties involved in a design are 
simulated by refined and accurate numerical 
models. In the end the use of “high fidelity” 
numerical simulation will also lead to overall 
cost reduction, since the analysis and design 
phase will be completed faster and only the 
reliability of the final configuration needs to be 
verified by structural testing. 

The light-weight shell structures used in 
aerospace applications are often buckling 
critical. The buckling load calculations are 
usually carried out by one of the many 
currently available finite element based 
computer codes [e.g., 3,4], In order to reduce 
computer execution time, buckling analyses 
are often done using only the small 
displacement stiffness matrix K 0 . This 

approach is used, despite the fact that the 
“initial stability problem” so formulated can only 
give physically meaningful answers if the 
elastic solutions based on K 0 (at least 
approximately) are identically equal to zero [5]. 

When the qualitative nature of the expected 
behavior is completely unknown, the stability of 
the structure must be investigated using the full 
tangent stiffness matrix Kj in order to 

guarantee accurate and reliable buckling load 
and buckling mode predictions. In order to 
discover the load level at which Kj ceases to 
be positive definite (that is, the load level when 
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buckling occurs), a step-by-step analysis 
procedure is needed. 

In addition, it is imperative (though often 
completely neglected), that at the beginning of 
any stability investigation, the accuracy of the 
discrete model used should be checked 
against available analytical or semi-analytical 
results. This step is part of a mandatory study 
needed in order to establish the dependence of 
the buckling load predictions on the mesh 
distribution used. Furthermore, as has been 
pointed out in the past by Byskov [6], if one 
carries out imperfection sensitivity 
investigations, which involve an extension of 
the solution into the postbuckling response 
region, further mesh refinement may be 
needed since the wavelength of the dominant 
large deformation pattern may often decrease 
significantly. 

Finally, whenever one is engaged in shell 
stability analysis it is especially important that 
one is aware of the possible detrimental effects 
of a whole series of factors, that have been 
investigated extensively in the late 1960s and 
the early 1970s. Thus for an accurate and 
reliable prediction of the critical buckling load 
of a real structure, one must account not only 
for the influence of initial imperfections [e.g., 
7,8] and of the boundary conditions [e.g. 9], but 
one must also consider the effects of stiffener 
and load eccentricity [e.g., 10] and the 

prebuckling deformations caused by the edge 
restraints [e.g. 1 1 ,12]. 

A recent test series of 5 composite shells 
carried out by Waters [1] at NASA Langley 
Research Center is used to illustrate how such 
a hierarchical approach to buckling load 
calculations can be carried out. The platform 
for the multi-level computations, needed for an 
accurate prediction of the critical buckling 
loads and a reliable estimation of their 
imperfection sensitivity, is provided by 
DISDECO [13]. With this open ended, 
hierarchical, interactive computer code the 
user can access from his workstation a 
succession of programs of increasing 
complexity. 


SOLUTION OF THE BUCKLING PROBLEM 


In the following it will be shown that with the 
help of DISDECO, the Delft Interactive Shell 
DEsign COde, the shell designer can study the 
buckling behavior of a specified shell, calculate 
its critical buckling load quite accurately and 
make a reliable prediction of the expected 
degree of imperfection sensitivity of the critical 
buckling load. The proposed procedure 
consists of a hierarchical approach, where the 
analyst proceeds step-by-step from the simpler 


(Level-1) methods used by the early 
investigators to the more sophisticated 
analytical and numerical (Level-2 and Level-3) 
methods used presently. 

Level-1 Perfect Shell Buckling Analysis 

The geometric and material properties of 
the 8-ply, composite shell with symmetrical lay- 
up of Ref. [1] are listed in Table 1 . 

Table 1. 

Geometric properties of NASA layered composite 
shell AW-CYL-1-1 [1] - [±45 / 0 / 90] s 


t total (= 

h) =0.039976 

in (=1.01539 

mm 

) 

L 

= 14.0 

in (=355.600 

mm 

) 

R 

= 7.99945 

in (=203.18603 

mm 

) 

Ell 

= 18.5111x10 6 

psi (= 12.7629x1 0 4 

N/mm 2 

) 

E 22 

= 1 .64x 1 0 6 

psi (= 1.13074x10 4 

N/mm 2 

) 

G 12 

= 0.8706x10 6 

psi (= 6.00257x1 0 3 

N/mm 2 

) 

v 12 

= 0.300235 





Note: Symmetrical lay-up with 8 plys of equal 
thicknesses (= 0.004997 in) 


Assuming a perfect shell (W = 0) and the 
following membrane prebuckling state 


W (o) =hW v =hAi2- 
c 

F (o),.Ehil 2 

cR 2 

where 


(D 


\ = 


N 


x . 


N C f 


and c 


a cf 

= ^/ 3(^- 2 


Eh M 

a C f - — — , N c ^ - a C f n 
cR 


V ) 


then the nonlinear equations governing the 
prebuckling state are identically satisfied and 
the linearized stability equations reduce to a 
set of equations with constant coefficients. It 
has been shown in Ref. [14] that by assuming 
an asymmetric bifurcation mode of the form 


W^ = h sin mu — cos — (y - x« x) (2) 

L R 

where 

m = k = number of axial half waves 
n = f. = number of circumferential full waves 
x« = Khot’s skewedness parameter 
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one can reduce the solution of the linearized 
stability equations to an algebraic eigenvalue 
problem. Notice that the eigenvalue X mnT 
depends besides on the wave numbers m and 
n also on Khot’s skewedness parameter , a 
real number. The critical load parameter X c is 
the lowest of all possible eigenvalues. Thus 
finding X c involves not only a search over the 

integer valued wave numbers m and n but the 
search has to be repeated over a range of 
possible positive and negative real numbers for 
x« . Using the Level-1 computational module 
AXBIF [14] a search over integer valued axial 
half-wave numbers m and over a range of 
possible positive and negative real numbers 
tK yielded the lowest eigenvalues listed in 

Table 2 for the specified circumferential wave 
numbers n. 

Table 2. 

Buckling loads of the NASA layered 
composite shell AW-CYL-1-1 
Buckling load map for the perfect shell 
using AXBIF [14] (N C(? = - 2238.325 Ib/in) 


n = 4 

A,™ =0.371613 

(m = 7, 

tk = 3.320) 

n = 5 

A™ = 0.370302 

(m = 7, 

IK = 2.680) 

CD 

II 

C 

A[P = 0.370029 

(m = 7, 

IK =2.178) 

n = 7 

A™ = 0.365992 

(m = 1, 

IK =0.011) 

n = 8 

A™ = 0.370131 

(m = 6, 

tk =1.681) 

n = 9 

0.371 980 

(m = 6, 

TK = 1488) 

n = 10 

A™ = 0.372880 

(m = 5, 

in 

CO 

T“ 

1! 

n = 11 

A™ = 0.375261 

(m = 4, 

T~ 

C\J 

II 

n = 12 

A™ = 0.369089 

(m = 1, 

tk = 0.579) 

n = 13 

X™ = 0.370076 

(m = 1, 

tk = 0.653) 

n = 14 

X™ = 0.376748 

(m = 1, 

tk = 0.734) 


Notice that besides the absolute minimum 
ofX?= 0.365992 at n = 7 there is a local 
minimum of Xc 1 = 0.369089 at n = 12. 

To facilitate the interpretation of the 
numerical results obtained, DISDECO provides 
the user with various graphical interfaces. 

Thus the results of the search for the critical 
(lowest) buckling load X c can be displayed in 

a contour map as shown in Fig. 1. Using 
membrane prebuckling the critical eigenvalue 
is (see also Table 2) 

= 0.365992 


with m=1 half-waves in the axial direction and 
n=7 full waves in the circumferential direction. 
In order to provide a quick overview of the 
distribution of eigenvalues, the values 
displayed in the contour plot are re-normalized. 
Thus in Fig. 1 the following re-normalized 
eigenvalues are plotted 

m _ ^mnx 
Pc 0.365992 

Notice that the critical buckling load can be 
calculated using a simple multiplication 

N? =*cNc e = 

0.365992 (-2238.325) = -819.209 Ib/in 

Level-2 Perfect Shell Buckling Analysis 

To investigate the effects of edge constraint 
and of different boundary conditions on the 
critical buckling load of the perfect shell 
(W = 0) one has to switch to the Level-2 

computational module ANILISA [15]. In this 
module the axisymmetric prebuckling state is 
represented by 

W< o) =hW v +hw 0 (x) 

(3) 

F<0) +r2, 0 (x)1 

It has been shown in Ref. [15] that with these 
assumptions the prebuckling problem is 
reduced to the solution of a single fourth order 
ordinary differential equation with constant 
coefficients, which always admits exponential 
solutions. Closed form solutions for simply 
supported and clamped boundary conditions 
have been published in the literature [16]. 

For anisotropic shells the linearized stability 
equations admit separable solutions of the 
form 

=h[w-|(x)cosn0 + W2(x)sinn0] 

(4) 

2 

pCO =^5t!_[f 1 ( x )cosn0 + f2(x)sinn0] 
c 

y 

where 0 = — . 

R 

Using a generalization of Stodola's method 
[17] first published by Cohen [18] the resulting 
nonlinear eigenvalue problem is reduced to a 
sequence of linearized eigenvalue problems. 
The resulting ordinary differential equations are 
solved numerically by a technique known as 
“parallel shooting over N-intervals" [19]. Notice 
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that by this approach the effect of edge 
restraint and the specific boundary conditions 
are satisfied rigorously. To find the critical load 
parameter A c an n-search must be carried 

out, whereby one must be careful to find not a 
local minimum but the absolute minimum. As 
can be seen from the results presented in 
Table 3 the n-search using membrane 
prebuckling and a rigorous satisfaction of SS-3 
(N x = v = w = M x =0) boundary conditions for 
the stability problem now yields a local 
minimum of A c =0.364715 at n = 7 and an 

absolute minimum of A c = 0.364370 at n = 1 1 . 
Table 3 

Buckling loads of the NASA layered composite 
shell AW-CYL-1-1 (N c ^ = -2238.325 Ib/in) 
Buckling load map for the perfect shell using 
ANILISA [1 5] (B.C. N x = v = w = M x = 0) 


Prebuckling: 

Membrane 

Nonlinear 

CD 

II 

C 

A™ = 0.371920 

A" 1 =0.347378 

n = 7 

A™ =0.364715 

A? =0.337088 

n = 8 

A™ = 0.372592 

A? =0.339701 

n = 9 

A™ =0.371460 

A^ 0.330957 

o 

n 

c 

A™ = 0.367479 

Ag'= 0.329163 

n = 11 

A™ = 0.364370 

A^ 0.328594 

n = 12 

A,™ = 0.364551 

Ac =0.330271 

CO 

n 

c 

A™ = 0.368528 

Ac =0.334044 

n = 14 

A™ = 0.376314 

Ac =0.339388 


The most accurate Level-2 solutions are 
obtained when one employs a rigorous 
nonlinear prebuckling analysis. As can be seen 
from the results listed in Table 3, for this 
particular shell the critical buckling loads with 
nonlinear prebuckling are always lower than 
the corresponding results obtained using a 
membrane prebuckling analysis. Specifically, 
the local minimum of A^ = 0.337088 at n = 7 
is about 8% lower, whereas the absolute 
minimum of A^ =0.328594 at n = 11 is about 
11% lower. Notice that the critical load N c can 
be calculated easily by multiplying the lowest 
eigenvalue A c by the normalizing factor 
No/ 1 = -2238.325 Ib/in yielding 

N c = A C N C/ = -735.500 Ib/in (n = 1 1) 


In fig. 2 the critical buckling modes using 
membrane and rigorous nonlinear prebuckling 
are depicted. Notice that the solutions with 
nonlinear prebuckling differ significantly from 
the ones obtained using membrane 
prebuckling, especially at n = 11 where one 
observes a typical edge buckling type 
behavior. 

Level-3 Perfect Shell Buckling Analysis 

To verify the earlier predictions the finite 
difference version [20] of the well known shell 
analysis code STAGS [21] will be used. Due to 
the slightly skewed buckling pattern predicted 
by the Level-1 and Level-2 computations one 
is forced to model the whole shell. 

Initially a convergence study must be 
carried out in order to establish the mesh size 
needed for accurate modeling of the buckling 
behavior of the shell in question. For this 
purpose the asymmetric bifurcation from a 
nonlinear prebuckling path option was used, 
whereby the earlier results obtained with the 
Level-2 module ANILISA listed in Table 3 
serve as a reference. 

In the convergence study, at first, for a fixed 
number of mesh points in the axial direction 
(NR = 161) the number of mesh points in the 
circumferential direction (NC) was increased 
until the bifurcation load approached a 
horizontal tangent. As can be seen from Fig. 3 
the results converge to a limiting value from 
below at about NC = 201. Next, for a fixed 
number of mesh points in the circumferential 
direction (NC = 201) the number of rows (NR) 
was varied. This time convergence is from 
above and as can be seen from Fig. 3 the 
horizontal tangent is reached at about NR = 
201 . 

Using a mesh of 161 rows and 201 columns 
(a model with 99268 D.O.F.'s and a maximum 
semi-bandwidth of 795) and SS-3 boundary 
conditions (N x =-N 0 ,v = W = M X = 0) the 
following 3 lowest eigenvalues were obtained: 

4 1 * = 0.327759 (n = 1 1) -> 

N^ = A^N^ = -733.630 Ib/in 

A^ = 0.328652 (n = 10) -> 

N< 2 > = A^Vl^ =-735.631 Ib/in 

4 3 > =0.329421 (n = 1 2) -* 

N(?> = 4 3 >N c , = -737.352 Ib/in 
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Details of the critical buckling mode are 
displayed in Fig. 4. Notice that the sequence of 
the 3 lowest buckling loads and the 
corresponding buckling modes agree closely 
with the predictions obtained with the Level-2 
module ANILISA (see also Table 3). 


IMPERFECTION SENSITIVITY STUDY 


That initial imperfections may decrease the 
load carrying capacity of thin-walled shell 
structures is by now widely known and 

accepted. However, in order to calculate the 
effect of initial imperfections one must know 
their shape and amplitude, an information that 
is rarely available. 

In the absence of initial imperfection 

measurements, as a first step one must 
establish whether a given shell-loading 

combination is imperfection sensitive, and if 
the answer is positive to estimate how 
damaging certain characteristic imperfection 

shapes are. 

Single Axisymmetric Imperfection 

Based on Koiter’s pioneering work on the 
effect of initial imperfections [7,22] the simplest 
imperfection model consists of a single 
axisymmetric imperfection 


where 


. x 


w o M = t r h^ i cos in - 

A, C j — A, L 


fo(x) = 


X (1 + B 2 ia 2 ) Eh 3 


, X 


^Cj 2apA22 C 


£l cosiit— 


( 7 ) 


-I (1 + B2ia 2 ) 2 

=-{afDl1+ 

* 2a, A22 


Notice that the linearized stability equations 
become now a set of equations with variable 
coefficients. The reduced wave number aj 

* 

and the normalized stiffness coefficients A 22 > 
* * 

B21 and D11 are all listed in Ref. [14]. 

It has been shown in Ref. [14] that by 
assuming an asymmetric bifurcation mode of 
the form 


=hsinmn— cos— (y - x^x) (8) 

L R K 

where 

m = k = number of axial half waves 
n = t = number of circumferential full waves 
= Khot’s skewedness parameter 


W=h^-| cosing (5) 

where i is an integer denoting the number of 
half-waves in the axial direction and ^ is the 

amplitude of the axisymmetric imperfection 
normalized by the shell wall-thickness h. 

If one assumes that both the axial load and the 
boundary conditions are independent of the 
circumferential coordinate, then the prebuck- 
ling solution will also be axisymmetric, a fact 
that simplifies the solution considerably. 

Level-1 Analysis of Axisymmetric 

Imperfection 

Neglecting the effect of the prebuckling 
boundary conditions the nonlinear equations 
governing the prebuckling state admit the 
following axisymmetric solutions 

Vy(°) = hW v + hw 0 (x) 

(6) 

F (o),_s£i Xy 2 +( (x) 

cR 2 y 0W 


a Galerkin type approximate solution yields for 
the eigenvalues (read, buckling loads) X of the 
problem a characteristic equation in the form of 
a cubic polynomial 

^ _ (^mnt +2^cj - Ci^i5j = 2m)^ (9) 

+ (2^mnx + ^Cj + (^2 - ^1 )^l$i=2m }^Cj ^ 

~{^mnx + ^2^iSi=2m + (^3 -648^)^ }X. 2 =0 
where 

Si=2m =1 ifi = 2m ; 

= 0 otherwise 

5j m = 1 if i = m 
= 0 otherwise 

and the constants C-|,C2 .••• ace listed in 
Appendix C of Ref. [14]. 
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Here it must be remembered that one will 
only get any noticeable degrading influence of 
the assumed axisymmetric imperfection if f-j 

is negative and if the coupling condition i=2m 
is satisfied. The physical explanation for this 
can be found in Koiter's 1963 paper [22]. 
Furthermore, in order to obtain the smallest 
real root of Eq. (9), for a given axisymmetric 

imperfection ^ an n-search must be carried 

out. It should also be noticed that the terms 
involving the Kronecker delta 5j = 2m are a H 

linear in f-j, and thus they dominate the 

buckling behavior of the shell with 
axisymmetric imperfection. 

Assuming that the most likely axisymmetric 
imperfection of the steel mandrel used to lay- 
up the NASA composite shell AW-CYL-1-1 is 
given by 


W = h^ 1 cos2jr^ (10) 

the Level-1 DISDECO computational module 
AXBIF generated the solid curve shown in Fig. 
5. Notice that the curve is re-normalized by 
A.™ = 0.366892 , the critical Level-1 buckling 
load of the perfect shell computed using AXBIF 
[14] with membrane prebuckling for =0.0 
and n=7. Notice also that an initial imperfection 
amplitude equal to the wall thickness of the 
shell (^-| = -1 .0) generates a “knockdown 
factor” of Pc = X c /^c* = 0.486 , resulting in the 
following rather low buckling load 

X c = 0.486^ = 0.486(0.366892) = 0.178310 

N c = ^c^c^ 

0.1 7831 0(-2238.325) = -399. 1 1 5 Ib/in 

with = 0.0 and n = 7 full waves in the 
circumferential direction. 

Level-2 Analysis of Axisymmetric 

Imperfection 

Since the external loading, the boundary 
conditions and the assumed initial imperfection 
are axisymmetric, therefore the prebuckling 
solution will also be axisymmetric. It has been 
shown in Ref. [23] that by assuming 

= hWv + hw 0 (x) 

( 11 ) 

F(0) =^- { -^y 2 + R2f o (x)} 


the solution of the nonlinear partial differential 
equations governing the prebuckling state can 
be reduced to the solution of a single fourth 
order ordinary differential equation with 
constant coefficients, which can be solved 
routinely. 

For anisotropic shells the resulting 
linearized stability equations admit separable 
solutions of the form 


W™ =h[w-)(x)cosn6 + W2(x)sinn6] 

(12) 

2 

p(1) = — — [f-|(x)cosne + f2(x)sinn0] 
c 

y 

where 0 = —. 

R 

Solution proceeds as outlined on Ref. [23]. 
Using an updated version of the Level-2 
computational module ANILISA [24] and SS-3 
(N x =v = w = M x =0) boundary conditions 

one obtains the results presented in Table 4. 
Notice that a rigorous nonlinear prebuckling 
analysis was used and an n-search was 
carried out for each specified axisymmetric 

imperfection amplitude ^ . 


The values of Table 4 are plotted as the 
dashed curve in Fig. 5. A comparison of the 
results obtained via the Level-1 module AXBIF 
(solid curve) and the Level-2 module ANILISA 
(dashed curve) shows that also in the case of 
axisymmetric imperfections a rigorous pre- 
buckling analysis should be used. Especially 
for very small initial imperfection amplitudes 

(|^|<0.1) the Level-1 predictions are 

inaccurate and overestimate the critical 
buckling load. Notice further that both curves 

have been normalized by X™ = 0.366892 , the 


critical Level-1 buckling load of the perfect 
shell computed using membrane prebuckling 
by AXBIF [14] for t k = 0.0 and n = 7. This way 

the effect of using a rigorous prebuckling 
analysis becomes easily discernible. 

It is interesting to see that for small enough 
initial imperfection amplitudes ( | < 0.07 , 

say) the critical buckling load of ihe shell is 
insensitive to the initial imperfection shape 
specified by Eq. (10). Notice that the critical 
buckling modes have n = 1 1 full waves in the 
circumferential direction, and as can be seen in 
Fig. 6, a somewhat skewed buckling mode 
shape which is dominated by edge buckling. 
However, larger initial imperfections 
(|^1 >0.07, say) force the shell to respond in 
another mode shape with n = 7 full waves in 
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the circumferential direction and with 
practically straight nodal lines. Interestingly 
enough now the critical buckling load of the 
shell is sensitive to the axisymmetric initial 
imperfection shape specified by Eq. (10) and 
an initial imperfection amplitude equal to the 
wall thickness of the shell (£,-] = -1.0) 
generates a “knockdown factor” of 
Pc =^cAc 1 =0.418. It predicts thus the 
following rather low buckling load 

X c = 0.41 8X™ = 0.418(0.366892) = 0.153361 
N c = ^c^c/? ” 

0.1 53361(-2238.325) = -343.271 Ib/in 


(remember STAGS defines W positive 
outward) the following critical bifurcation load 
was found 

N c = -333.449 Ib/in 

As can be seen from Fig. 7 the critical buckling 
mode has n = 7 full waves in the circum- 
ferential direction and no visible skewedness. 
The nondimensional bifurcation load of the 
shell with axisymmetric imperfection is for 

fi=to 

Np -333.449 

X r = — — = = 0.148973 

u Nc^ -2238.325 


with a very slight skewedness of the buckling 
pattern and n = 7 full waves in the 
circumferential direction. 

Table 4 

Buckling loads of the NASA layered composite 
shell AW-CYL-1-1 (N c<? = -2238.325 Ib/in) 
Axisymmetric imperfection using ANILISA [24] 
(B.C. : N x = v = w = M x = 0) 


The re-normalized bifurcation load is 
Pc = where N c/? = Xj? 1 N c/ > 

N c < 

thus 


nl _ X c _ 0.148973 

Pc ~ ~ 0.327759 

v c 


0.4545 


*1 

CO 

$i 

i nl 
A c 

0. 

0.328594 (n=1 1) 

-0.2 

0.299830 (n=7) 

- 0.01 

0.328671 (n=1 1) 

-0.3 

0.274773 (n=7) 

- 0.02 

0.328745 (n=1 1) 

-0.4 

0.251275 (n=7) 

- 0.05 

0.328939 (n=11) 

-0.5 

0.229976 (n=7) 

-0.06 

0.328994 (n=11) 

-0.6 

0.210876 (n=7) 

-0.07 

0.328332 (n=7) 

-0.7 

0.193827 (n=7) 

-0.08 

0.326652 (n=7) 

-0.8 

0.178655 (n=7) 

-0.09 

0.324860 (n=7) 

- 0.9 

0.165187 (n=7) 

-0.10 

0.322960 (n=7) 

- 1.0 

0.153263 (n=7) 


Level-3 Analysis of Axisymmetric 

Imperfection 

Recalling that since both the axial load and 
the boundary conditions are independent of the 
circumferential coordinate, therefore the 
prebuckling solution will also be axisymmetric, 
one can use once again the asymmetric 
bifurcation from a nonlinear prebuckling path 
option. By modeling the full shell the code can 
choose itself the critical number of full waves In 
the circumferential direction. No n-search must 
be carried out. Using a uniformly spaced mesh 
of 161 rows and 201 columns and the user 
written subroutine option WIMP to introduce 
the following axisymmetric imperfection 

W = 1.0hcos27t^ 


Notice that the Level-3 re-normalization is 
done using = 0.327759 , the critical 
Level-3 buckling load of the perfect shell 
computed using STAGS-A with nonlinear 
prebuckling and a 161x201 mesh. 

Notice also that the Level-2 ANILISA prediction 
(pc’ =0.4664, n = 7) agrees closely with the 
Level-3 STAGS-A prediction (p{J = 0.4545, 
n = 7) . The slight difference is partly due to 
the fact that ANILISA uses the Donnell type 
nonlinear shell equations, whereas STAGS-A 
employs the higher order Marlowe-Fliigge 
equations. 

Single Asymmetric Imperfection 

The effect of a single asymmetric initial 
imperfection can be investigated either by 
solving the full nonlinear response problem or 
by employing the well known Lyapunov- 
Schmidt-Koiter [7] reduction technique. When 
investigating the degrading effect of a single 
mode asymmetric imperfection 

W = tf 2 sinmn-^cosn^- (13) 

where m and n are integers denoting the 
number of axial half-waves and the number of 
circumferential full waves, respectively, 
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instability occurs at the limit point of the 
prebuckling state in the generalized load- 
deformation space. Assuming that the 
eigenvalue problem for the critical (lowest) 
buckling load A c will yield a unique asym- 
metric buckling mode W^\ then for an 
imperfect shell (f 2 * 0) the shape of the 
generalized load-deflection curve in the vicinity 
of the bifurcation point A = A c is given by the 
following asymptotic expansion 

(A- A c )^ = A c a^ 2 + A c b^ 3 + ... 

_ (14) 

-A c 0^2 -(A-A c )^ 2 +0(^2) 

Expressions for the postbuckling coefficients 
“a” and “b" and the imperfection forms factors 
“a” and “P” are derived in References [25,26], 
If the limit point is close to the bifurcation point, 
then the maximum load A s that the structure 
can carry prior to buckling can be evaluated 
from Eq. (14) by maximizing A with respect to 
E, . For cases where the first postbuckling 
coefficient “a” is zero, this analysis yields the 
modified Koiter formula [26] 

(1-p s ) 3/2 =|V-3a 2 b[1-i(1- Ps )p 2 |(15) 
where p s = A s / A c . 

Notice that, if the second postbuckling 
coefficient “b" is positive, Eq. (15) has no real 
solutions. Thus the buckling load of the 
specified shell-loading combination is not 
sensitive to small asymmetric initial 
imperfections of the shape given by Eq. (13). 
If, however, the second postbuckling 
coefficient “b” is negative, the equilibrium load 
A decreases following buckling and the 
buckling load of the real structure A s is 
sensitive to the asymmetric initial imperfection 
specified by Eq. (13). 

Level-1 Analysis of Asymmetric 

Imperfection 

For the composite shell under investigation, 
as can be seen from the partial results listed in 
Table 2, there are many eigenvalues only 
slightly higher than the critical one of 
X c = 0.365992 for m = 1 , n = 7 and 
= 0.01 1 . Hence, strictly speaking, the 
proposed form of the perturbation expansion 
given by Eqs. (14) is not applicable, since the 
nonlinear interaction between the many nearly 
simultaneous eigenmodes is not accounted for. 
Thus the following results, where one 
considers the eigenfunctions corresponding to 


certain critical eigenvalues chosen one at the 
time, can at best give an indication as to the 
severity of the expected imperfection 
sensitivity. 

Assuming initially an asymmetric 
imperfection affine to the critical buckling mode 
of the perfect NASA composite shell AW-CYL- 
1-1 as computed by the Level-1 computational 
module AXBIF (see also Table 2) 

W = hf2Sin7c— cos — (y- 0.01 lx) (16) 

L R 

and using the Level-1 computational module 
BFACT to carry out the initial postbuckling 
analysis yields the following results 

X c = V? = 0.365992 (m = 1, n = 7, x K = 0.01 1) 

b = - 0.048844 a = P = 1.0 

Substituting these values into Eq. (15), one 
can plot the degrading effect of an asymmetric 
imperfection of the shape given by Eq. (16) as 
a function of its amplitude ^2 • As can be seen 
from Fig. 8 an initial imperfection amplitude 
equal to the wall thickness of the shell 
(^2 = 10) generates a “knockdown factor” of 

p s = Ag/A™ =0.541 , resulting in the following 
rather low buckling load 

X s = 0.541 A™ =0.541 (0.365992) = 0.1 98002 

Ng = A s N c /> = 

0.1 98002 (-2238.325) = -443.192 Ib/in 

Notice that the imperfection form factors 
"a” and 11 P" are identical equal to 1.0 because 
BFACT uses membrane prebuckling to 
calculate the necessary first and second order 
fields and the assumed asymmetric 
imperfection shape of Eq. (16) is affine to the 
buckling mode. Please notice that in Fig. 8 the 
collapse load is re-normalized by 

A™ = 0.365992 , the critical Level-1 buckling 
load of the perfect shell computed using 
membrane prebuckling by AXBIF [14] for 
= 0.01 1 and n = 7 . 

Level-2 Analysis of Asymmetric 

Imperfection 

To investigate the effects of edge-constraint 
and/or different boundary conditions on the 
imperfection sensitivity of the critical buckling 
load of the NASA composite shell AW-CYL-1-1 
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one has to switch to the Level-2 module 
ANILISA [15] and run its postbuckling analysis 
option. In this module, as described earlier, the 
axisymmetric prebuckling state is represented 
by Eqs. (3), the buckling modes by Eqs. (4) 
and the postbuckling state by 

W^' = h[w a (x) + wp (x) cos n0 + Wy (x) sin n0] 

(18) 

2 

p(2) =^— _ [f a ( x ) + fp (x) cos 2n0 + fy (x) sin 2n0] 

where 0 = y / R . Details of the computational 

procedures used are reported in Refs. [15,23]. 

Next, let us assume that the specified 
asymmetric imperfection is affine to the critical 
buckling mode obtained by the rigorous Level- 
2 perfect shell analysis discussed earlier 

W = h^ 2 [w-|(x)cos11 0 + W2(x)sin11 0] 

where 0 = y/R and the component functions 
w-|(x) and W 2 (x) are shown in Fig. 2b. 
Running ANILISA with rigorous prebuckling 
and SS-3 boundary conditions (N x = -N 0 , 
v = w = M x = 0) yields the following results 

X c =0.328594 (n = 11) 

b = -0.37605 ; a = 0.46663 ; p = -0.22174 

Using Eq. (15) to plot the degrading effect 
of the asymmetric imperfection specified by 
Eq. (18) as a function of its amplitude ^2 one 

obtains the results displayed in Fig. 9 as a 
solid line. Obviously the fact that for an 
imperfection shape affine to (similar to) the 
buckling mode with an amplitude of ^2 =1-0 

one obtains a negative load carrying capacity 
is unrealistic. 

Here one must remember that Koiter’s 
Sensitivity Theory is asymptotically exact, that 
is, it yields accurate predictions for sufficiently 
small imperfections, whereby what is 
sufficiently small may vary from case to case. 
Also, Eq. (15) was obtained by using the 
perturbation expansion given by Eq. (14), 
where terms of order (J£) are neglected. As 

can be seen from the dotted curve plotted in 
Fig. 9, by using more advanced computational 
modules such as COLLAPSE [27], where a full 
nonlinear solution is used and terms up to and 
-2 

including order {Zfe ) are kept, one obtains 
more reasonable predictions. 


Notice that up to about ^2 = 0 - 3 the 

asymptotic predictions from ANILISA and the 
nonlinear results of COLLAPSE agree very 
closely. Thus one can say that in this case the 
range of validity of the asymptotic solution is 

0 > > 0.3 . 

DISCUSSION OF THE RESULTS 

When comparing and analyzing the results 
obtained sofar it is important to keep in mind 
that all Level-1 and Level-2 solutions are 
based on approximate representations of the 
unknown functions. As pointed out in the 
previous sections Level-1 solutions use a 
single term double Fourier series 
approximation to reduce the solution of the 
stability problem, formulated in terms of partial 
differential equations, to algebraic eigenvalue 
problems. The effect of edge restraint is 
neglected (one uses a membrane prebuckling 
solution) and the assumed field functions 
satisfy approximately SS-3 (N x =-N 0 . 
v = w = M x =0) boundary conditions. 

Level-2 solutions eliminate the y- 
dependence by a truncated Fourier 
decomposition in the circumferential direction. 
The resulting system of nonlinear ordinary 
differential equations are solved numerically, 
whereby both the specified boundary 
conditions and the effect of edge restraint are 
rigorously satisfied. Thus by this approach the 
only approximation is that one represents the 
variation of the solution in the circumferential 
direction by a single harmonic with n full 
waves, whereby an n-search is used to 
establish which wave number is the critical 
one. The Level-2 module ANILISA can also be 
used to investigate the effect of using different 
boundary conditions. In Table 5, the results for 
four different boundary conditions are 
presented. As expected the fully clamped C4 
boundary conditions has the highest critical 
buckling load. The increase in load carrying 
capacity with respect to the weaker SS3 
boundary conditions is about the same as for 
an isotropic shell of similar characteristic 
dimensions (same L/R and R/t ratios) of Ref. 
[28], 

The Level-3 solutions are based either on a 
2-dimensional finite difference or finite element 
formulation. In both cases, if one uses the 
appropriate meshes, one can obtain rigorous 
solutions where all nonlinear effects are 
properly accounted for. The only real problem 
with Level-3 type solutions is that for each 
problem one must establish the appropriate 
mesh size. Coarse meshes yield inaccurate 
solutions. What is coarse depends on the 
particular problem under investigation. Thus, 


9 

American Institute of Aeronautics and Astronautics 



for a general nonlinear solution a convergence 
study must always be carried out. 

Using a hierarchical simulation platform 
such as DISDECO (Delft Interactive Shell 
DEsign Code), where the analyst has at his 
disposal computational modules of different 
level of sophistication, such a convergence 
study can be carried out relatively quickly and 
accurately. In Table 5 a summary of the results 
obtained in this study is presented using 
normalized variables. In Table 6 the same 
results are repeated but this time the imperfect 
buckling loads are printed as re-normalized 
variables p . Looking at the first column, where 
the critical buckling loads and the critical 
buckling mode shapes of the perfect shell are 
listed, one sees that using the Level-3 code 
STAGS-A one must indeed use a relatively fine 
mesh (161 rows and 201 columns) in order to 
obtain an accurate prediction of the critical 
buckling load. Remember, all preceding 
computational modules are based on Donnell 
type anisotropic shell equations, however, 
STAGS-A uses the more accurate Marlowe- 
Flugge type equations. Usually the use of a 
more refined theory implies a lower buckling 
load. The value of the normalization factor 
used, N c /=Eh 2 /cR, is printed in the 
heading of the table. The lowest critical 
buckling load was found using STAGS-A and 
modeling the whole shell with a mesh 
consisting of 161 rows and 201 columns. The 
asymmetric bifurcation from a nonlinear 
prebuckling path option with SS-3 boundary 
conditions yielded 

N" 1 =0.327759 (-2238.325) 

= -733.630 Ib/in (n = 1 1) 

Considering now the effect of different 
types of imperfection shapes a second 
normalization is introduced, whereby the new 
normalization factors are chosen such that for 
vanishingly small imperfections the normalized 
variable p approaches unity (i.e. 1.000). The 
only exception to this rule is the case of the 
axisymmetric imperfection 

W = -h^-j cos2n-jj 

where by using as the normalization factor 

= 0.366892 (the asymmetric bifurcation 
perfect shell buckling load with membrane 
prebuckling, = 0.0 and n = 7) in the limit as 
£l 0 , p c approaches the value 


p c = = Q ' 3285 — = 0.89561 5 - 0.896 

KC i 0.366892 

A C 

This value represents the effect of edge 
restraint. See also the results of Table 4 and 
Fig. 5. Notice that using membrane 
prebuckling, thus neglecting the effect of edge 
restraint, the normalized buckling load of the 
perfect shell is 1 .000. 

Turning now to the effect of asymmetric 
imperfections, from the results listed in Table 5 
it is evident that for an imperfection amplitude 
equal to one wall-thickness the range of 
validity of the asymptotic solutions is 
exceeded, and the predictions of Koiter’s 
imperfection sensitivity theory computed by 
ANILISA [15] are no longer valid (see also Fig. 
9). On the other hand, there seems to be good 
agreement between the results obtained by 
COLLAPSE [27], a Level-2 computational 
module which computes a nonlinear solution 
based on a two modes approximation, and the 
STAGS-A [20] solution obtained sofar for the 
two modal imperfections and the affine 
imperfection listed in Table 5. In general the 
more accurate STAGS-A, Level-3 solutions are 
slightly lower than the Level-2 COLLAPSE 
solutions, with exception of the asymmetric 
imperfection affine to the perfect shell buckling 
mode. The reason for this anomaly lies in the 
fact that this imperfection triggers more than 
one circumferential harmonic close to the limit 
point (as can be seen in Fig. 66 of Ref. 29). 


CONCLUSIONS 

By relying on a series of theoretical results of 
various degree of sophistication published in 
the literature, the hierarchical approach used in 
this paper has resulted in a series of buckling 
load predictions of increasing accuracy. It was 
shown that in order to be able to arrive at a 
reliable prediction of the critical buckling load 
and to make an estimate of its imperfection 
sensitivity which can be used with confidence, 
one must proceed step by step from simple to 
more complex models and solution proce- 
dures. 

In particular one can state, that in order to 
predict the critical buckling load accurately and 
to make a reliable estimate of its imperfection 
sensitivity, the nonlinear effects caused by the 
edge restraint conditions must be included in 
the analysis. Any solution procedure which 
fails to account for these effects, should be 
suspect of having provided incorrect results. 

The most approximate of the here 
described analyses, the Level-1 solutions 
which neglect the effects caused by the edge 
restraints, can still be used to great advantage 
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to establish the approximate behavior of a shell 
subjected to the specified external loading. 
However, depending on the value of the 
prebuckling stiffness, resulting from the 
different types of wall constructions used, the 
solutions may be either conservative or 
nonconservative. 

As can be seen from the results shown in 
Tables 5 and 6, the buckling load of the 
composite shell AW-CYL-1-1 is sensitive to all 
the initial imperfection shapes investigated. For 
a more specific prediction of the final collapse 
load, the final goal of a “High Fidelity Buckling 
Load Analysis”, one would have to carry out a 
refined Level-3 analysis including measured 
values of all the significant generalized 
imperfections such as the traditional shell-wall 
imperfections, variations in the shell-end or 
loading surface geometry and especially for 
composite shells variations in the shell-wall 
thickness distribution. It has been shown in 
Ref. 30 that such an approach yields very good 
agreement between the predicted collapse 
load and the experimental buckling load. Such 
a paper is in preparation. 
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number of axial half waves 


Contour map of the buckling loads for the NASA composite shell AW CYL 1 1 

rhoc — No / Ncl-bor where N c I — bor*l ombdoc *Ncl“ — B 1 9 209 lb/m 



contour 

level s 
1 1 .00 
2 1.10 
3 1 .20 
4- 1 .50 
5 2.00 


number of circumferential full waves n 


Fig. 1 Distribution of buckling loads based on Level-1 membrane prebuckling 
analysis - NASA composite shell AW-CYL-1-1 
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buckling, n=7 
5 




c. Nonlinear pr 


X" 1 =0.329 







a. General view of the critical buckling mode 


NASA composite shell AW-CVL- 1 - 1 
: — trace critical buckling mode — Nc = 


— 161 rows x 201 columns — 55 — 3 B.C. 

-733.640 Ib/in — delta -teta = 360 degrees 



b. Axial trace of the critical buckling mode 

Fig. 4 Buckling mode of the axially compressed layered composite shell AW-CYL-1-1 
SS-3 B.C. -N x = -N 0 ,v = W = M x = 0; N C f = -2238.325 Ib/in; STAGS-A results 
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Fig. 5 Imperfection sensitivity for axisymmetric imperfection under axial compression 
(SS-3B.C.: N x = -N 0 ,v = W = M x = 0 ; =-2238.325 Ib/in) 


o. Prebuckling Shape 


b. Buckling Component Modes | 

= -0.07, 0.328332 (n =7) 

W - X 

Fig. 6 Critical buckling modes for axisymmetric imperfections - — = cos2tc— 



b. Buckling Component Modes 


^1 = -0.06, Ju- 0.328994(n=1 1) 
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NASA composite shell AW-CYL-1 -1 — 161 rows x 20 1 columns — SS 3 B.C 

critical buckling mode at Nc = -333.4-A9 Ib/in - Wba r/h = 1 . 0 - c os < 2 ♦ p I * x/L) 



Fig. 7 General view of the critical buckling mode - axisymmetric imperfection 


imperfection sensitivity of NASA composite sHelt 01 SS3 B.C 
Wbcjr = h * x i fc>a r2 * s in ( pi * x/L) * c os f "7 * (y 0.01 1 



normolizod imperfection amplitude — xtbor2 


Fig. 8 


Level-1 imperfection sensitivity calculation using BFACT [14] 
_ x 7 

Asymptotic theory: W = h^sinrc— cos— (y- 0.01 lx) 

L R 
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Table 5 SUMMARY OF RESULTS - Axially compressed NASA layered composite shell AW-CYL-1 -1 
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Table 6 SUMMARY OF RESULTS - Axially compressed NASA layered composite shell AW-CYL-1-1 

t'W = -2238.325 Ib/in - laminate f±45/0/90L 
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Boundary conditions: SS3: 




